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Abstract 

This article proposes a model of bone remodeling that encompasses mechanical and 
electrical stimuli. The remodeling formulation proposed by Weinans and collaborators 
was used as the basis of this research, with a literature review allowing a constitutive 
model evaluating the permittivity of bone tissue to be developed. This allowed the 
mass distribution that depends on mechanical and electrical stimuli to be obtained. 
The remaining constants were established through numerical experimentation. The 
results demonstrate that mass distribution is altered under electrical stimulation, 
generally resulting in a greater deposition of mass. In addition, the frequency of 
application of an electric field can affect the distribution of mass; at a lower frequency 
there is more mass in the domain. These numerical experiments open up discussion 
concerning the importance of the electric field in the remodeling process and propose 
the quantification of their effects. 
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Background 

Bones provide mechanical stability to the human body and are a source of minerals for 
metabolism [1]. Bones have been studied extensively from the mechanical and mineral 
standpoint, and in terms of functionality [1,2]. From a mechanical point of view, bones 
can be adapted to loads on trajectories of stress through mineral apposition, which is 
due to the action of osteoblasts [1-4]. Furthermore, they reabsorb minerals when the 
mechanical stimulus is sufficiently low, as it is unnecessary to maintain structure [2]. 
Reabsorption is directed by osteoclasts. Osteoblasts and osteoclasts are the primary 
cells involved during bone remodeling that are stimulated by the action of mechanical 
strain sensors, for example, osteocytes [2]. These three cell types play an important 
role during the processes of replacement, maintenance and modeling of bones [1]. 

Following the work of Meyer during the nineteenth century, Wolff [5] proposed a the- 
ory of trabecular bone architecture, which assumes that trajectories of high mechanical 
stress form the trabecular bone. In 1987, Frost [6-8] suggested an adaptive mechanism of 
bone mass as a function of mechanical stress. Consequently, several bone remodeling 
algorithms have been developed including those proposed by Frost [8], Pauwels [9], Kum- 
mer [10], Cowin [11-13] and Hegedus [14], which predict the formation of bone structure 
from internal mechanical loads studied in terms of stress and strain. 

From mechanical models of bone remodeling, sophisticated studies have been car- 
ried out concerning the processes of apposition and reabsorption during bone 
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turnover, and particularly concerning the distribution of mass in the femur [15,16], hip 
replacement [17,18], and dental implants [19]. Generally, these studies have been phe- 
nomenological. Therefore, researchers have made significant efforts to include mathem- 
atical models and the role of cell biology and biochemistry in the remodeling process, 
resulting in research that begins at the microscopic level, concerning the effects of 
basic cellular remodeling units (BMU, Basic multicellular units) during tissue replace- 
ment [20,21]. From the perspective of BMU, important work was initiated at the bio- 
chemical and mechanical level concerning the effects of cracks [22], cell cycles 
throughout adult life [23], active molecules within each cell [24] and the spatial distri- 
bution of each BMU [25]. With these important advances in the understanding of bone 
remodeling, researchers in the field increasingly turned to the study of other biophys- 
ical stimuli that can affect this process. Most models have not taken account of the 
physical-chemical phenomena of tissue mechano-transduction. For this reason, new 
investigations that allow the piezoelectric and electrokinetic behavior of the bone to be 
studied were undertaken [26] . 

A clinical study demonstrated that a local electromagnetic field accelerates the heal- 
ing process after bone fracture [26]. Therefore, an article by Demiray and Dost [27] 
began new research concerning the effect of the electromagnetic field on interior injury 
to bone. In another article, Ramtani [26] presented a mathematical model relating to 
the benefit of the electric field in the reparair and maintenance of the solid matrix of 
bone. Furthermore, there are studies concerning the electrical behavior of bone tissue 
during the production of electric fields, and external electrical flow. Fukada and Yasuda 
[28] demonstrated that bone exhibits piezoelectric behavior, i.e. mechanical stress cre- 
ates electric polarization (the indirect effect) and an external electric field causes strain 
(the converter effect). In addition, the properties of bones that produce piezoelectric 
potentials have been determined [29-33]. These data led to the development of math- 
ematical models that include the effect of electromagnetic fields during the repair 
[34,35] and remodeling [36] of bone. For example, Qu and Yu [34] developed a math- 
ematical model (no spatial dimension) of the remodeling process and healing under the 
effect of mechanical loads and the use of electric charges. In this model, the higher the 
voltage applied to a bone after fracture, the lower the percentage of bone damage and 
micro damage in the few days after the stimulus. Similarly, during osteoporosis an elec- 
tric field increases bone density over time. Huang et al. [37] established a hypothesis 
concerning the biological and biochemical pathways that activate cells, particularly 
osteocytes, during the imposition of an electric field. Furthermore, Qu and Yu [38] pro- 
posed a mathematical model that included mechanical loads and electromagnetic 
effects during the process of bone remodeling. 

To date, there have been no phenomenological models concerning bone remodeling that 
have been tested and compared with purely mechanical models. Therefore, this article pro- 
poses a new electro-mechanical model relating to bone remodeling. To test its performance, 
various numerical experiments were carried out and compared with previous mechanical 
models. The electric model constants were obtained from relevant literature and numerical 
experimentation. From these assumptions it was concluded that the electric field can affect 
the distribution of mass, which originates from the remodeling process, under mechanical 
effects only. Using the remodeling model of Nackenhorst [39] as a starting point, it was 
demonstrated that the electric field can increase bone density and accelerate the process of 
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apposition. Therefore, the model proposed herein can be used as a basis for further work 
concerning electrical effects in the maintenance of bone. 

Methods 

The electro-mechanical model 

The electro-mechanical model of bone remodeling that involves mechanical and elec- 
trical stimuli can be written, hypothetically, as follows (1): 

~ = gmechip, W mec (p)) +gelect(p, W e l ect (c(p,f))) (1) 

at 

Whereg mec h(p, W mec (p)) is the well known mechanical stimulus described by Weinans 
[39], which depends on tissue density (p(x, y, z, t)), and the work carried out by mech- 
anical stress (W mec (p)) and g e kct{p, W e i ect (e(p,f))) is the electrical stimulus that depends 
on density, frequency and the work carried out by the electric field ( W e i ect (c(p))). Dur- 
ing this first approach, we consider that the two stimuli are added to determine the 
bone remodeling process. However, we will develop each of the terms that determine 
the electro-mechanical model throughout this article. 

The mechanical model 

Following the remodeling mechanical model described in [4,39], the density variation 
over time depends on the mechanical stimulus that exists at every spatial point of the 
bone, which can be written as in [39] (2): 

gmechip, W mec (p)) = h - 1^ (2) 

Where pis the bone density at each point in space (p(x,y,z, t)), W(p)is the energy- 
strain per unit of volume due to mechanical stress, k\ is a constant and W re j m is the 
strain-energy (per unit of volume) of reference that sets the threshold for which they 
will perform deposition (W(p)/W re f m > 1) or absorption (W(p)/W re f m < 1) of the tissue 
[39], in the presence of mechanical stress. It should be noted that the energy-strain 
depends on the density and is given by (3): 

W(p) = ~e T C(p)£ (3) 
2p 

Where e is the strain, in Voigt notation, of the strain tensor given in (4): 

e=[£n £22 y 12 } T ( 4 ) 
That is a function of the displacements given by (5): 
u = [ui u 2 ] T (5) 
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In addition, C(p) is the matrix of linear elasticity. The matrix C(p) contains the Pois- 
son module, which is usually considered constant, and Young's modulus that depends 
on the density by expression [4] (7): 

E{p)=Ap" (7) 
where A is a constant and n establishes a relationship of power density that has been 
uncovered through experimentation [39]. 

By manipulating equation (7) we can obtain a dimensionless form of the density. This 
is easier to use with the aim of determining the modulus of elasticity. Multiplying the 
right side of (7) by (pq/Pq) produces (8): 

E(p) = Apl (£-) "=Eo(£Y = E 0 A" (8) 



Where Eq = Ap^ and A = p/p 0 are the elasticity modulus and the dimensionless density 
ratio, respectively. Therefore, the linear elasticity matrix can be expressed as (9): 

C(p) = A"Q> (9) 

Where Co is the matrix of linear elasticity with constant coefficients, which depends 
onEo and v only, and is given in the case of plane stress by: 
1 v 0 

(10) 



Co 



1 



v 1 0 

o o i(i-v) 



Thus, the strain energy per unit of volume (3) can be expressed as (11): 

W(p) = ^AVQe = - = - U mec (11) 

2p p V 2 J p 

Where U mec is the strain energy at each instant of time, which is calculated with the 
initial constant of the remodeling [39] problem only. Replacing these equations in (2), 
and with some algebraic manipulation, we produce (12): 

^■K-^- 1 ) (12) 

where we can define U re f m = p Q W re f m . Therefore, we produce the following equation for 
the density ratio (13): 

&«*=*i(A"- 1 ^-l > ) (13) 



U r 



The momentum equation that establishes the internal stresses of a body is given by 
[40] (14): 

V«a + b = 0 (14) 
Where the stress is given by (15): 
a = C(p)e = A"Coe (15) 
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Electrical model 

This article proposes the inclusion of a hypothetical electrical term that can determine, 
in part, the process of bone remodeling. Therefore, the contribution of this stimulus 
can be written as (16): 

geleAp, W elect (c(pJ))) = k2 Wd ^ Pj)) (16) 

Where € {p,f) is the electrical permittivity of bone tissue and depends on the dens- 
ity, (p) where frequency (/*), W e i ect (<E (p,f)) is the electrical energy per unit of volume, 
I<2 is a constant and W re / e is the electric energy (per unit volume) of reference. 

It should be noted that the term electric energy depends on permittivity, which in 
turn depends on the density and frequency. This energy term is given by (17): 

W elect (e (p,f)) =^e (p,f)E(x,y,z, if (17) 

Where ~E(x,y, z, t) is the electric intensity (electric field). However, e(p,f), electrical 
permittivity, depends on density and frequency, and is given by (18): 

e ( P J) = e 0 e r ( P J) (18) 

Where Go is the permittivity in the free space and G r (p,/) is the relative permittivity, 
which in turn is given by (19): 

e r (p,f)=B P m + S(f) (19) 

where B is a constant and m establishes a relationship of power density that can be 
proved through experimentation, and will be developed in the following sections. Fur- 
thermore, 8(f) is a function of frequency at which the electric field is applied. 

As with the mechanical case, we can manipulate equation (19) so that it is 
expressed in terms of relative density, multiplying the first term by to produce 
ipl/pt) (20): 

e r ( P ,f) = Bp™X m + 8(f) = e p X m + 6(f) (20) 

where <E p is a constant for the power model of relative permittivity. Therefore, sub- 
stituting (20) in (17), we obtain (21): 



W elect (e (pJ)) = ^e 0 {e p X m + S(f))-E(x, y ,z.,t) 2 

z PPo 

1 , . m ,2 1 



-e 0 (e p X m + 8(f))E{x,y,z,tf = — (e p X m + S(f))U dect 



(21) 



2 Po X p 0 X 
where U e iect = ^o^(x,y,z, if. Substituting equation (21) in (16) produces (22): 

( w r , MVl . + . {e p X m + S(f)) U dect 

gelectip, W e l ect (e (/>))) = h ~ ™ = *2 " , ^77 (22) 

We have chosen U re f e = p 0 W re f c as the reference value of electric energy. 
However, Gauss's law for an electric field, with no internal loads, is given by (23): 

V • D = p e (23) 



where D =e E and p e is the density of the electric energy. In turn, the electric field 
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can be expressed in terms of a quantity called electric potential or voltage, given 
by (24): 

E = -V0 (24) 
where0 is the electric potential. In summary, equation (1) can be written as (25): 

dp = k ( yn—\ H^_j\ + h M" + 8 M (25) 
dt \ Uref m J X U re f e 

Again, using the dimensionless form, we have (26): 

d\ _ ki A„_ t U mec \ | k 2 (g^ m + 8(f)) Uelect ^ 

dt p 0 \ U refm J Pq X U refe 

dX ( \n-\ Umec A , (gpA" + 8(f)) U elect 

~T — kmec I A — 1 I + Kelect ^ 7} 



where k mec = A:i//3 0 and k e i ec t = kijp^ are mechanical and electrical constants that define 
the conversion rate of bone remodeling, dependent on the mechanical stress and the 
electrical potential, respectively. 

Solution by the finite element method 

To solve equation (14) we take its differential form to its weak [40] form, to obtain 
(27): 

8eX n C 0 zdn - J 5u • bdO - J 8u • tdr = 0 (27) 
a n r 

Where 8e and 8u are weighting functions, and t = a • n is the traction on the 
boundary r that serves as a frontier to the domain D (see Figure 1). 

To calculate the approximate solution by finite element discretization, the displace- 
ment field approximates through [40] (28): 

w(x, t) = N r (x)a(t) (28) 

Where N r (x) is a row vector containing the shape functions used to approximate the 
displacement a(t) [40] in the nodes. Using the Galerkin method, we approximate the 
weighting functions in the same way as the displacement field, producing at the elem- 
entary level equation [40] (29): 

(J BX"CoB T dn e ^a+ I NbdQ e + B.C. = 0 (29) 

n' 

where B is the derivative operator (discrete version of operator Bi of equation (6)) that 
converts the displacements into strains (see [40]). This system of equations is com- 
pleted by applying the Neumann and Dirichlet conditions suitable for solving the elastic 
problem (see Figure 1). 

Similarly, for the electric case we can use the finite element method. For this 
purpose we use the weak expression of equation (23) in terms of electrical 
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elect 

Dirichlet 



r=r 




Neumann 



Neumann 



"p mec 
■•■Dirichlet 



Figure 1 Domain and boundary conditions. The acronym on the index indicates whether it is a 
condition for the electric (elect) or mechanical (mec) case. Note that there are two types of boundaries (not 
necessarily equal) for each equation, the electric field and displacement (for the mechanical stress). 



potential to produce (30): 

(Vw)(e p x m + S(j)){v<p)dn 



wp e dfl 



/ ((e p x m + S{f)){v<p)*n)dr 



r 
0 



(30) 



where w is the weighting function and n is the external normal to the domain of 
the problem. Again, for calculation of the approximate solution, the scalar field of 
the electric potential approaches through [40] (31): 

?(x,£) = N r (x)v(t) (31) 
Where v is the potential at each node of the element, and N(x) are the shape func- 
tions for scalar [40] problems. Weighting functions are chosen in the same manner as 
shape functions, which is through the method of (Bubnov-) Galerkin standard. There- 
fore, at the elementary level we have the following equation (32): 



B e (e p A m + S(f))Bjdn e 



N T p e d£2 + B.C. = 0 



(32) 



where 





dNi - 


dx 


dy 


dN 2 


dN 2 


dx 


dy 


dN nod 


dN nod 


. dx 


dy . 



(33) 



where nod is the number of nodes of each element. 
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Solution to the equation of relative density 

This article solves the equation of dimensionless density (or relative) by integrating 
equation (26) by the method of Euler. For this purpose we define (34): 



where At = t k+ \ — t k is the integration time interval and k refers to the evaluation of 
the variable A at a specific time, i.e.: = A(x, t k ). This method has been used exten- 
sively in the prediction of bone density through remodeling [42]. The forward Euler 
method is of the first order and has the disadvantage of being unstable for large time 
intervals [41]. 

Euler's method was implemented in FORTRAN and was coupled with the elastic and 
electrical problems. For implementation we used an approach based on element (with 
an elemental average) [39,43]. 

Computational model 

The computational example to be solved is a 0.1x0.1 m square plate with non-uniform 
load at the top and restrictions of movement vertically at the bottom as presented in 
Figure 2. This example has been extensively studied in numerical tests concerning bone 
remodeling algorithms [42-44]. For these values we built a mesh of 80x80 elements 
(see Figure 3) in the horizontal and vertical directions, respectively, and integrated 
these with a time step of At = 0.100. The total simulation time was 200. The initial 
condition was A = 1.0 and the limits used in the algorithm of bone remodeling were: 
0.0125<A<2.175. 

The constants for the mechanical energy 

For the mechanical case we used an initial Young's modulus E Q , a Poisson's modulus v 
and an initial dimensionless density Ao with values 64 MPa, 0.3 and 1.00, respectively 
[44]. The dimensionless parameters of the density equation were obtained from Nack- 
enhorst (1997) [39,44] and were n = 2, k mec = 03Y15days~ l and U refm = 800Pa. 

The mesh was produced using bilinear quadrilateral elements and four points of 
Gauss integration [40]. 

The constants for electric energy 

In the electric case, there are numerous articles [45-48] that determine the major elec- 
trical properties of bone. In particular, graphs of the electric permittivity in function of 
frequency and density are provided in [46], as presented in Figures 4 and 5. Figure 4 
presents data concerning bone density as a function of the relative permittivity on the 
frequency of 50 KHz. The article [46] used 40 bone tissue samples, and used the meas- 
urement of the density as a function of the permittivity. From the data, the mathemat- 
ical model proposed for permittivity as a function of tissue density is given by a power 




(34) 



Therefore, the forward Euler method is defined as [41] (35): 



(35) 
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equation, expressed as (36): 



e w (p, 50KHz)=Bp" 



(36) 



Where m = 1.5486 and £ = 1050.0; this is a function of the density of bone tissue. 
With the aim of using dimensionless density, we multiplied and divided by (oq/Pq) to 
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Relative permittivity 

Figure 4 Bone density as a function of the relative permittivity at a frequency of 50 KHz. Taken from 
[46]. 



produce (37): 

€ W (p) = 1050.0^ 5486 A L5486 = GpA 1 ' 5486 ( 37) 

Where Gp = 1050. OpJ 5486 and depends on the initial density that considers the com- 
putational model. 

Figure 5 presents the relative permittivity in terms of the function of the frequency 
for the density of the femoral head of p = 03737g.cm~ 3 . With a greater frequency, per- 
mittivity decreases. To complement equation (37) (33) we must add a term that allows 
us to calculate the difference between the relative permittivity at a frequency of 50 kHz 
and any other frequencies used in the model; this is (38): 

<%)(/■) = e r{f) (03737g.cm-\f) - e r{50Khz] (03737g.cm-\f = 5QKhz) 

S(f)if) = e r(f) (03737g.cm- 3 J) - 228.681 {M) 

where & r (50Khz) (0- 3737g. cm~ 3 ,f = 50Khz) = 228.681 . For its part, the function 




10 I 

100 1k 10k 100k 1M 10M 

Frequency (Hz) 



Figure 5 Graphic representation of relative permittivity (dimensionless) as a function of frequency 
(in Hz) for various densities of bone tissue. FC: Femoral head; FMC: femoral medial condyle, FLC: femoral 
lateral condyle FTM: femoral greater trochanter. Taken from [46]. 
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€ r (fl (0.3737g\ cm 3 ,f) can be obtained from Figure 5, so we get: 

log w [e r(f) {0.3737g.cm- 3 J)] = 0.1407 [lo glo (/0] 2 - 1.8936 log 10 (f) 

+ 8.1505 (39) 

In summary, the equation for the electric permittivity depending on the density and 
frequency is given by (see Figure 6): 

e r (p,f) = e^A 1 ' 5486 + e r(J] - 228.681 (40) 

For the remaining constant changes we made variations of k e i ec and used U re f e = 
800-Pa, as in the mechanical case. 



Results 

Figure 7 presents various examples that were run for a frequency f = 50 KHz, and with 
variation in the constant K e iect- For these examples, a potential difference at various 
contour lines was imposed, according to the domain that was established in Figure 2. 
In the first example, a voltage of 100 was placed on the right side and no voltage at the 
bottom (Figure 7a). In the second example, a voltage of 100 V was placed on the right 
side and no voltage on the left (Figure 7b). In the third and fourth examples, the volt- 
age was placed on the top and bottom, respectively (Figures 7c and 7d). In the other 
contours we placed null Neumann conditions. 



Relative permittivity in function of Density and frequency 
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In the first example, the mechanical condition demonstrated in Figure 2 coupled with 
an electric field as observed in Figure 7a were used. The results are presented in Figure 8. 

Figure 8 presents the results for the relative density value of \ using an approach 
based on element, with an average at an elementary level of the mechanical stimulus, 
electrical and density values. In all figures (Figures 8, 9, 10, 11 and 12) black represents 
the maximum density value (2,175) and white represents the lowest value (0.0125). It 




Figure 8 Results for the first example: the electric boundary condition in Figure 7a. The figure 
demonstrates the results for various values of k e | ect . 
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should be noted that with At = 0.1 the pattern obtained is similar for all values of k e i ect . 
It is observed that near the area of imposition of the mechanical stress there is a similar 
density distribution to a chess board, and far from the loading area there is the forma- 
tion of three well-defined columns reminiscent of the cortical bone (Figure 8a and 8b). 
In cases where the constant increases (k e i ect = 1.4xl0 7 and above) there is the appear- 
ance of a high density area in the lower right region, from where a new column starts 
in the top domain (Figs 8c, d, e and f). In this area, near the lower right structure, there 
is the distribution and space formation of a chess board. In addition, in the upper part 
there are empty areas that generate ramifications in each of the columns supporting 
the load. Note that the imposition of the electric field defines a new topological struc- 
ture in the domain, as presented on the right-hand side of the simulation results. This 
new structure is an additional column, which is generated by the potential and a sup- 
port area of higher density in the lower right corner. 

Figure 9 presents the results for various values of K e iect- hi this case, a voltage of 
100 V was imposed on the right side and a null voltage on the left side. Null Neumann 
conditions were imposed on remaining contours (see Figure 7b). It is noted that there 
were changes in topology; columns 2 and 3, from left to right in Figures 9c and 9d, are 
thicker and closer to one another. In addition, upper branches of greater density were 
created and an additional branch that begins from the last column to the right. In the 
case of Figure 9e, the formation of a non-defined region of the "chessboard" type on 
the right side of the domain is presented. 

Figures 10 and 11 present results for various values of k e i ect . In these cases a voltage 
of 100 V was imposed at the top and bottom, respectively. Null Neumann conditions 



Garzon-Alvarado et al. Theoretical Biology and Medical Modelling 2012, 9:14 
http://www.tbiomed.eom/content/9/1 /1 4 



Page 15 of 17 




were imposed on the remaining contours. Note that in these cases the mass distribu- 
tion increases as k e i ect increases. In the case of Figure lOe, the density increases in the 
upper part of the domain, so that the chessboard becomes more continuous in the cen- 
tral part of the domain. In Figures lOe and lie the columns that are formed in each 
simulation have higher bandwidths than previous cases. 

In the second example there is a variation of the frequency, while the value of k e i ect 
remains constant at 7.0xl0 7 . With the voltage configuration of Figure 7d, this is with a 
lower voltage of 100 V. Note in this case that at low frequencies the density and the 
amount of tissue deposited are greater than at high values. Figure 9a demonstrates that the 
frequency generates a total deposit of tissue. In figure 12b the formation of a high density 
topology and continuing formation is apparent. Figures 12c, 12d and 12e present results 
comparable to those observed in previous cases. However, in figure 12c, the columns are 
wider than in figures 12d and 12e. 

Discussion and conclusions 

In this article several numerical examples were developed concerning bone remodeling 
of the plate during mechanical stress, assuming the imposition of an electric field in 
the domain. To calculate the mechanical and electrical stimulus of remodeling, and the 
evolution of density, the elemental approach was utilized. This pioneering article 
includes the electrical effect, previously designed by Weinans et al. [4], in a model of 
bone remodeling. 

Comparable with previous articles, the results of the study presented herein dem- 
onstrate in the mechanical case formation of the columnar zone (of high density) 
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in the area remote from the load and formation of the trabecular zone (of low 
density) in the area close to the load [17]. The results are similar to those obtained 
by Weinans et al. [4], Fernandez et al. [49] and Chen et al. [17]. When applying 
an electric field there is an increase in bone density and an alteration in the top- 
ology of the distribution of mass in the domain. In general, there is greater bone 
mass apposition in the domain. Therefore, the columns developed by the mechan- 
ical stress increase in size due to the electric field. In addition, a greater number 
of columns and localized compact zones may be observed. In the formulation the 
effect of electrical frequency has been included to allow increased apposition of 
mass at low frequency to be observed. 
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